#!/usr/bin/env python

"""
StirlingModel.py: An ODE model for oxygen uptake kinetics
"""

from numpy import trapz

class StirlingModel:
	# Model parameters variables
	A = float(0.0)
	B = float(0.0)
	C = float(0.0)
	E = float(0.0)
	VdotO2min = float(0.0)
	VdotO2max = float(0.0)

	# Calculated data
	Time_Grid = []
	Demand = []
	VdotO2aerobic = []
	VdotO2anaerobic = []
	VO2aerobic = float(0.0)
	VO2anerobic = float(0.0)

	def dVdotO2dt( self, D, VdotO2 ):
		# Solve the Stirling model's equation at one particular set
		# of conditions.
		dVdt = self.A * ( VdotO2 - self.VdotO2min )**self.B * ( self.VdotO2max - VdotO2 )**self.C * ( D - VdotO2 )**self.E 
		return dVdt

	# Solve function.  If you provide a demand curve, it will take this
	# model and solve it for the course of the demand curve.
	def solve( self, t, D, VdotO20 ):
		self.Time_Grid = t
		self.Demand = D
		self.VdotO2aerobic = self.solveAerobic( VdotO20 )
		self.VdotO2anaerobic = self.solveAnaerobic()
		self.VO2aerobic = trapz( self.VdotO2aerobic, x=self.Time_Grid )
		self.VO2anaerobic = trapz( self.VdotO2anaerobic, x=self.Time_Grid )

	def solveAerobic( self, VdotO20 ):
		nt = len( self.Demand )
		VdotO2 = []
		VdotO2.append( VdotO20 )

		for n in range(0, nt-1):
			dt = self.Time_Grid[n+1] - self.Time_Grid[n]
			w1 = dt * self.dVdotO2dt( self.Demand[n], VdotO2[n] )
			w2 = dt * self.dVdotO2dt( 0.5*(self.Demand[n]+self.Demand[n+1]), VdotO2[n]+0.5*w1 )
			w3 = dt * self.dVdotO2dt( 0.5*(self.Demand[n]+self.Demand[n+1]), VdotO2[n]+0.5*w2 )
			w4 = dt * self.dVdotO2dt( self.Demand[n+1], VdotO2[n]+w3 )
			VdotO2.append( VdotO2[n] + (1.0/6.0) * ( w1 + 2.0*w2 + 2.0*w3 + w4 ) )
		return VdotO2

	def solveAnaerobic( self ):
		nt = len( self.Demand )
		VdotO2anaerobic = []

		for n in range(0, nt):
			if ( self.Demand[n] > self.VdotO2min ):
				VdotO2anaerobic.append( self.Demand[n] - self.VdotO2aerobic[n] )
			else:
				VdotO2anaerobic.append( 0.0 )
		return VdotO2anaerobic

	def getDemand( self ):
		return self.Demand

	def getAerobic( self ):
		return self.VdotO2aerobic

	def getAnaerobic( self ):
		return self.VdotO2anaerobic
	
	def getAerobicTotal( self ):
		return self.VO2aerobic

	def getAnaerobicTotal( self ):
		return self.VO2anaerobic

	# Constructor
	def __init__( self, A, B, C, E, VdotO2min, VdotO2max ):
		# This creates a new model equation without a demand curve.  This
		# model can be thought of as a person, and the demand curve is
		# just a workout or race they did.
		self.A = A
		self.B = B
		self.C = C
		self.E = E
		self.VdotO2min = VdotO2min
		self.VdotO2max = VdotO2max
